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Impurity in a Bose-Einstein condensate: study of the attractive 
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We investigate the properties of an impurity immersed in a dilute Bose gas at zero temperature 
using quantum Monte-Carlo methods. The interactions between bosons are modeled by a hard 
sphere potential with scattering length a, whereas the interactions between the impurity and the 
bosons are modeled by a short-range, square-well potential where both the sign and the strength 
of the scattering length b can be varied by adjusting the well depth. We characterize the attractive 
and the repulsive polaron branch by calculating the binding energy and the effective mass of the 
impurity. Furthermore, we investigate structural properties of the bath, such as the impurity-boson 
contact parameter and the change of the density profile around the impurity. At the unitary limit 
of the impurity-boson interaction, we find that the effective mass of the impurity remains smaller 
than twice its bare mass, while the binding energy scales with jm, where n is the density of 

the bath and m is the common mass of the impurity and the bosons in the bath. The implications 
for the phase diagram of binary Bose-Bose mixtures at small concentrations are also discussed. 

PACS numbers: 


I. INTRODUCTION 

The polaron problem is a paradigmatic topic in con¬ 
densed matter physics: it concerns the effect of the quan¬ 
tum fluctuations of the surrounding medium on the prop¬ 
erties of an impurity immersed in a bath. The first for¬ 
mulation of the problem is due to Landau and Pekar [T] 
in their study of the motion of electrons in polar crys¬ 
tals. By using a variational approach valid in the limit 
of strong coupling it was shown that the electron be¬ 
comes eventually trapped in the potential created by the 
self-induced deformation of the lattice. Frohlich [2] pro¬ 
posed an effective Hamiltonian, which describes the cou¬ 
pling between charged impurities and longitudinal opti¬ 
cal phonons of the lattice, providing the standard model 
description of the polaron problem. The ground-state en¬ 
ergy of the Frohlich Hamiltonian was calculated by Feyn¬ 
man [3] within a variational approach based on path in¬ 
tegrals which yields an upper bound that nicely interpo¬ 
lates between the weak-coupling perturbative result and 
the strong-coupling Landau-Pekar prediction. After a 
large amount of theoretical work [4] spanning more than 
four decades, an exact solution of the model by means 
of diagrammatic Monte Carlo methods was finally pre¬ 
sented in Refs. 1311]. Remarkably, both the results ob¬ 
tained for the polaron ground-state energy and the ef¬ 
fective mass are in very good agreement with Feynman’s 
findings. 

Important generalizations of the Frohlich Hamiltonian 
include the Holstein model on a lattice [7| , electrons cou¬ 
pled to acoustical phonons in a crystal [8| and, more 
recently, impurity particles immersed in a dilute Bose- 
Einstein condensed (BEC) gas [3HT3]. 

In the context of ultracold atoms the polaron concept 
has received large attention in its fermionic version, i.e. 
an impurity coupled to a Fermi sea m- Thanks to the 


use of a Feshbach resonance the s-wave scattering length 
between the impurity and the fermions of the bath can be 
tuned at will and experiments have probed various prop¬ 
erties of attractive and repulsive Fermi polarons both in 
3D [T3IfT3 and in 2D ranging from the weak to 

the strong-coupling regime and including the polaron- 
molecule transition. 

Bose polarons, which involve a bath consisting of a 
BEC and therefore are directly related to the original 
Frohlich model, have also been realized and their dy¬ 
namics experimentally investigated [T5II20) . However, so 
far, there have been no studies exploiting Feshbach res¬ 
onances to increase the strength of inter-species inter¬ 
actions, nor measuring basic polaron properties such as 
their binding energy, lifetime and effective mass. On the 
theoretical side, the self-localization of Bose polarons was 
investigated using mean-field approaches [Mil 1211111] 
as well as Feynman’s variational method applied to the 
effective Hamiltonian describing the impurity mm- 
Starting from the Frohlich Hamiltonian other studies 
have focused on the calculation of the radio frequency 
response of the polaron [53], and of its binding energy 
and effective mass using renormalization group |25j and 
diagrammatic Monte Carlo [53] methods. A more micro¬ 
scopic approach based on the T-matrix approximation 
was used in Ref. [27] where various quasiparticle proper¬ 
ties are calculated for both attractive and repulsive Bose 
polarons close to a Feshbach resonance. Similar results 
are also obtained in Ref. [55] by means of a variational 
ansatz for the wave function of the bath-impurity sys¬ 
tem. Finally, three-body correlations were explicitly in¬ 
cluded in the theoretical treatment both at the level of 
perturbation theory |29] as well as within a variational 
approach [3D]i giving rise to a significant lowering of the 
binding energy of attractive polarons. 

In this paper we address the problem of Bose polarons 
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using a fully microscopic, non perturbative approach, 
consisting in the quantum Monte Carlo (QMC) method. 
This numerical technique can provide exact results for 
the ground-state energy and the effective mass of the 
impurity as a function of the parameters of the Hamil¬ 
tonian describing the inter-species and intra-species in¬ 
teraction potentials and the density of the bosonic bath. 
We model these interactions using a hard-sphere poten¬ 
tial for the inter-boson repulsion and both a purely repul¬ 
sive hard-sphere and an attractive square-well potential 
for the impurity-boson interaction. In particular, the lat¬ 
ter model allows one to investigate situations where the 
impurity-boson s-wave scattering length is either posi¬ 
tive or negative giving rise to the ground-state attrac¬ 
tive and excited-state repulsive branches of the polaron. 
Our analysis is limited to the case where the mass of the 
impurity is equal to the one of the Bose particles in the 
medium, but generalizations to include different mass ra¬ 
tios can be easily implemented within the same method. 

We investigate the properties of the Bose polaron both 
along the attractive and the repulsive branch. We find 
that for small values of the ratio \b\/a of the impurity- 
boson to the boson-boson scattering length our results 
for the binding energy and the effective mass are in good 
agreement with second-order perturbation theory based 
on a Frolich-like Hamiltonian describing the coupling be¬ 
tween the impurity and the bath. At the unitary point of 
resonant impurity-boson scattering (h = ±oo) the bind¬ 
ing energy is found to scale with where n is 

the density of the bath and m is the common mass of the 
particles and the impurity. We notice that this behavior 
is similar to the Fermi polaron case where the binding en¬ 
ergy is proportional to the Fermi energy of the bath m- 
The effective mass ratio m*/m ranges from values close 
to one in the weak-coupling limit up to values which re¬ 
main smaller than two close to the resonant point. We 
find no evidence of the self-localization of the polaron, 
which in studies based on the Frohlich model is signaled 
by an abrupt increase of the effective mass as the cou¬ 
pling strength exceeds a critical value. We believe that 
this wrong prediction has to do with the inadequacy of 
the effective Frohlich Hamiltonian in the description of 
the pairing mechanism which takes place close to the res¬ 
onance where the impurity and one boson from the bath 
can form a bound state. 

We analyze the structural properties of the bosonic 
bath by calculating the contact parameter which char¬ 
acterizes the short-range behavior of the impurity-boson 
pair correlation function. The knowledge of how particles 
in the bath are distributed around the impurity enables 
us to evaluate the distortion of the density profile pro¬ 
duced by the impurity. Within the attractive square-well 
model, we find a pronounced peak in the density close to 
the impurity both on the attractive and on the repulsive 
branch. This peak is a result of the pairing induced by 
the impurity-boson potential. It is a short-range feature 


that can not be accounted for by the Frohlich Hamilto¬ 
nian which can only describe long-range distortions of 
the density profile. 

An important point to analyze is related to the exis¬ 
tence of few-body bound states in vacuum, such as three- 
body Efimov-like states and deeper bound states with 
more than three particles. At the resonant point we cal¬ 
culate the energy of the deepest bound state with three 
and more particles (i.e. the impurity plus two or more 
bosons) finding evidence that such state exists only up 
to six particles {i.e. the impurity plus five bosons). Re¬ 
markably, the energy of these self-bound states is in abso¬ 
lute value much smaller than the polaron binding energy 
which involves the contribution from a large number of 
particles in the bath. One should notice that the values 
we obtained for the ground-state energy of the cluster 
states as well as the size of the largest cluster greatly 
depend on the details of the hard-sphere boson-boson 
potential used in the simulations. However, we believe 
that the results for the polaron binding energy at unitar- 
ity are universal and only depend on the gas parameter 
na^ of the bath and the mass ratio between the impurity 
and the bosons. 

The structure of the paper is as follows. In Sec. H we 
first address the single-polaron problem by introducing 
the model Hamiltonian (subsection H-A), and by review¬ 
ing the perturbation treatment leading to the Frohlich- 
type Hamiltonian. Here, we also derive the results for 
the polaron binding energy and effective mass valid in 
the weak-coupling limit (subsection H-B). Finally, in sub¬ 
section H-C, we briefly review the DMC method and we 
discuss the different trial wave functions used to describe 
the attractive and repulsive polaron branch. The results 
along the two branches concerning binding energy, effec¬ 
tive mass, density profiles and contact parameter are pre¬ 
sented in subsection H-D. Furthermore, subsection H-E 
contains a discussion of these results specific of the reso¬ 
nant point for the impurity-boson scattering. In Sec. HI 
we report on calculations of the binding energy of few- 
body states in vacuum at the unitary point and on the 
side of the resonance where a two-body bound state ex¬ 
ists. In Sec. IV, we generalize the problem to many im¬ 
purities obeying Bose statistics and we use DMC simula¬ 
tions to validate the perturbative equation of state in the 
limit of small concentrations bearing some consequences 
for the phase diagram of binary mixtures. Conclusions 
are finally drawn in Sec. V. 


II. SINGLE IMPURITY 

A. Model Hamiltonian 

We consider a system of one impurity immersed in a 
dilute gas of N Bose particles at T = 0 described by the 
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following Hamiltonian 


.2 N 


H = -- 


2m B 


2mj 


i—l i<j 

N 


( 1 ) 


Here, the first two terms represent the kinetic and 
the interaction energy of the bosonic bath consisting 
of particles of mass ms and interacting through the 
two-body potential Vb, which depends on the distance 
rij = |ri — r^l between a pair of bosons. Furthermore, 
— h'^/{2mi)'\7^ is the kinetic energy of the impurity with 
mass TO/ denoted by the coordinate vector Tq, and V7 is 
the boson-impurity potential depending on the distance 
fia = Itq — f/l between the impurity and the /-th particle 
of the bath. The inter-boson potential Vb is modeled by 
the hard-sphere (HS) interaction 


Vsir) 


-boo r < a 
0 r > a , 


( 2 ) 


where the diameter a coincides with the s-wave scatter¬ 
ing length. The impurity-boson interaction, instead, is 
modeled by either a purely repulsive hard-sphere poten¬ 
tial 


where = 2 to_r(Vo — |e6|)/?i^. We also notice that the 
value KqRq = 7r/2 corresponds to the unitary limit of the 
impurity-boson interaction where the scattering length 
b diverges and the binding energy e/ vanishes. In the 
present study we consider values of the range Rq that are 
small compared to the interboson scattering length with 
a ranging between 5 to 20i?o. We expect that for such a 
short-range potential the value b of its scattering length 
is the only relevant parameter for all polaron properties. 

We restrict the analysis of the Hamiltonian Q to the 
case where the impurity and the bosons in the bath have 
the same mass: to/ = to^ = 2m r = to. The strength 
of the inter-boson interactions is determined by the gas 
parameter na^ involving the bosonic density n = N/V, 
whereas the intensity of the impurity-boson coupling is 
given in terms of the ratio \b\/a of the two scattering 
lengths. 


B. Perturbation theory 

The problem of a single mobile impurity in a Bose gas 
can be thoroughly investigated using perturbation the¬ 
ory at least in the weak-coupling regime. The approach 
is based on the treatment of the bath within the Bogoli- 
ubov approximation of a dilute Bose gas described by the 
Hamiltonian 


Vi{r) = VPir) = 


- 1-00 r < b 


^0 r > b , 
or an attractive square-well (SW) potential 

^-Vq r < Ro 


Vi{r) = V7^(r) = 


0 


r > Rn . 


( 3 ) 


( 4 ) 


The latter is characterized by a range Rq and a depth Vq 
(Vb > 0) chosen such as to yield the value b of the scat¬ 
tering length. This is determined by the transcendental 
equation 


tan{KoRo) 

KqRq 


( 5 ) 


where = 2TO/jVb/?i^ in terms of the reduced mass 
TO// = mjmB / (mj -b to^). In the case of the repulsive 
potential V/^ the s-wave scattering length is always pos¬ 
itive, whereas for the square-well potential in Eq. ([^ the 
value of b can be either positive or negative depending 
on KqRq. In particular, we consider values in the range 
0 < KqRq < TT, Corresponding to either no bound state 
{KqRq < •7r/2) or one bound state {KqRq > n/2) in the 
two-body sector. In this latter case the molecular binding 
energy e/ is obtained from the equation 

tan(Kj?o) _ h 

kRq i?o\/2TO//|eb| ’ 


Hb = Eb + ^ etalak . (7) 

k 

Here, Eb is the ground-state energy of the bosonic par¬ 
ticles 

= • <*> 

where g = ° is the inter-boson coupling constant. 

The operators a^, are the annihilation/creation op¬ 
erators of quasiparticles related to the bosonic particle 
operators Uk, aj^. through the standard transformations 

Oik = U/;Uk U/gut.k 

aj(. = Ukol - UfcO-k , (9) 

with coefficients u1 = 1 + and UkVk = 

— 1^. The elementary excitation energies are given by 
the Bogoliubov spectrum 

Efc = y(e°)^ + ^guQel , (10) 

where is the dispersion of free particles and Uq 

is the density of condensed particles. 

At the mean-field level the interaction energy between 
the impurity, located in position Tq,, and the bath is de¬ 
scribed by the following expression 

Hint = gsi J dr n{r)S{r - Ta) , (H) 
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involving the density of bosons and the interspecies cou¬ 
pling constant ggj = ^ proportional to the boson- 

impurity s-wave scattering length b. In momentum space 
the above interaction Hamiltonian can be recast in the 
form 


Hint 


9Bin (12) 

V ^ q#0 


where the Bogoliubov approximation a|j.ak+q — 
t/noV{uq + Vq){aq + OfLq) of the density fluctuations 
in terms of quasiparticle operators 0 has been used. 
By applying perturbation theory to the Hamiltonian 
Hb + Hint, as it was done for example in Ref. [3T] for 
a Bose-Fermi mixture, one finds the following result for 
the ground-state energy of the system of N bosons plus 
one impurity 

Eo = Eb + gn (-+ , (13) 

\a 3y7r J 


valid to order in the boson-impurity coupling 
strength m EH- To the same order b^ one can also 
calculate the effective mass of the impurity obtaining the 
result [3 [33] 


= 1-H 


64 

45-y7r 


(14) 


ratio, m*/m ^ 1, which arises from the cloud of bosonic 
quasiparticles dressing the impurity. We would like to 
stress, however, that the Frohlich Hamiltonian HB+Hint, 
where the bath and interaction term are given respec¬ 
tively by Eq. 0 and ( [T^ ), is an effective low-energy re¬ 
duction of the microscopic Hamiltonian Q. Whether it 
captures the relevant physics of an impurity immersed 
in a Bose condensate when the coupling to the bath is 
strong is questionable and should be analyzed with care. 
This problem will be investigated in the remaining part of 
the article using quantum Monte Carlo methods that are 
particularly suitable to treat strongly correlated systems 
in a non perturbative manner. 


C. Quantum Monte Carlo method 


We use the diffusion Monte Carlo (DMC) method 
which aims to solve the many-body Schrodinger equa¬ 
tion in imaginary time t = it/Ti for the distribution func¬ 
tion /(R, r) = ■!/'T(R)d>(R, r), where 4)(R, r) is the wave 
function of the system of N bosons plus one impurity and 
'0 t(R) is a trial wave function of the particle coordinates 
R = (ri,..., r^v, fa) used for importance sampling. The 
time-dependent Schrodinger equation can be written as 


9/(R,t) 

dr 


-i4V|/(R, r) + D\/n ■ [F(R)/(R, r)] 
[EL{-R)-Eo]fi-R,T) , (15) 


Apart from the trivial first-order contribution gBin to 
the ground-state energy, the second-order corrections in 
Eqs. (13l-(14) scale in terms of the same dimensionless 
parameter t/ruE^, which should be much smaller than 
unity to ensure the validity of the perturbation approach. 
We also notice that the perturbative corrections to the 
energy and the effective mass in Eqs. (13)-(14| diverge 
when the inter-boson scattering length a tends to zero. 
This feature indicates the instability of the ideal Bose gas 
towards clusterization around the impurity and points 
out the crucial role played by the repulsive interaction 
between the bosons. 

The Hamiltonian (121 has the general form of the 
Frohlich polaron Hamiltonian describing an impurity 
coupled to a bath of non interacting bosonic quasiparti¬ 
cles P] . This analogy was first exploited in Ref. [5| where 
the intriguing problem of the “self-localization” of the po¬ 
laron was addressed in a fashion similar to the Landau- 
Pekar description of electrons in ionic crystals [T|. In 
Ref. [13 the Jensen-Feynman variational scheme is ap¬ 
plied to both the strong and the weak-coupling regime 
of the effective Frohlich Hamiltonian. Both these stud¬ 
ies predict that for > 0.7 the impurity “self- 

localizes” inside the potential well produced by its dis¬ 
tortion of the bosonic density. In terms of the effec¬ 
tive mass m* of the impurity the self-localization phe¬ 
nomenon corresponds to a large enhancement of the mass 


in terms of the so called local energy Ei(R) = 
(i7^/!T(R)/V’T(R) and quantum drift force F(R) = 
(2 Vr'!/'t(R))/'0t(R)- In the above equation D = 
/{2m) plays the role of a diffusion constant and Eq is 
a reference energy. The formal solution of the equation 
is given by 

/(R', t + 5t)= J dR G(R', R, (5t)/(R, r) , (16) 


where one introduces the Green’s function 
G(R',R, (5t) = (R'|e“‘^'”^|R) describing the time 

evolution governed by the Langevin operator 
C = —-I- DVr • F -I- {El — Eq). If the short-time 
dependence of G(R',R, Jr) is known for sufficiently 
small 6t, the asymptotic solution for large times, 
/(R,t —>■ oo), can be obtained by iterating Eq. ( [l^ ) for 
a large number of time steps 5t. On general grounds, the 
initial time distribution /(R, 0) = '(/'t(R) X)™ Cn‘bn(R) 
is expanded in terms of the eigenfunctions $„(R) of the 
system corresponding to the eigenenergies E„. For a 
system of bosons, provided the coefficient Cq does not 
vanish, the solution at large times of Eq. (151 is given 
by the expression 


/(R, T oo) = coV’t(R)$o(R) , (17) 

and is proportional to the nodeless ground state 4>o(R) 
with energy Eq. The value of the ground-state energy is 
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determined from the condition of keeping the probability 
distribution /(R, r) stationary at large times or, more 
conveniently, from the average of the local energy E£(R) 


J dR /(R, T — oo)El{R) 
J dR /(R, T —)■ oo) 


(18) 


Apart from statistical errors, the DMC method allows 
one to calculate the exact ground-state energy of a sys¬ 
tem of Bose particles. Importantly, the energy estimate 
obtained using the DMC technique with importance sam¬ 
pling is to a large extent independent of the detailed 
shape of the trial wave function as long as '0'r(R-) is posi¬ 
tive definite. If the function ■07’(R) changes sign in some 
regions of the configuration space, the DMC algorithm 
guided by this trial function yields, instead of the ground 
state, the lowest-energy eigenstate compatible with the 
nodal constraint fixed by 07 " 

The general form of the trial wave function used in the 
present study is given by 


N 

0t(R) = fiina) n feinj), (19) 

i=l i<j 


where the functions fs and // describe, respectively, 
inter-boson and impurity-boson two-body correlations. 
The functional form of the inter-boson term is con¬ 
structed from the two-body scattering solution of the 
hard-sphere potential in Eq. ([^ 


1. Repulsive branch (& > 0): fj is constructed from the 
zero-energy scattering solution of the potential @ , or¬ 
thogonal to the bound state existing at 6 > 0 for two 
particles: 

r r < Ro 

= Ro<r<R (22) 

R<r<L/2. 


Here R is a matching point and for r > R the function 
// goes to a constant reached at L/2 where fj{r = 
L/2) = 0. The coefficients A, B and C ensure the 
continuity of // and of its first derivative at the points 
Rq and R. 


2. Attractive branch {b < 0): // is constructed in the 
same way as for the repulsive branch of Eq. (22), with 
the only difference that the scattering length b is neg¬ 
ative. 


3. Attractive branch (6 > 0): We use the solution of the 
two-body bound state with energy eb given by Eq. ([^ 


Ro<r<R (23) 
R<r<L/2, 

where k} = m{Vo — |e{,|)/?i^. The coefficients A, B 
and C again ensure the continuity of fj and of its 
first derivative at the potential range Rq and at the 
matching point R. 


fsir) = 


0 

sin[fe(r—a)] 


r < a 

a < r < L/2 , 


( 20 ) 


where the value of the wave vector k is chosen such that 
the first derivative of the function vanishes at half of 
the size L = of the cubic simulation box: f'^{r = 
L/2) = 0. This condition ensures that the Jastrow factor 
(201 is compatible with the periodic boundary conditions 


used in the simulation. 

For the impurity-boson correlation function // we use 
instead different forms depending on the type of poten¬ 
tial Vi, hard sphere or square well, and on the polaron 
branch, repulsive or attractive. 

• Hard-sphere potential: We use the two-body scattering 
solution of the potential (|^ 


In all the above three cases, the values of the matching 
point R and of the parameter a are optimized by mini¬ 
mizing the variational energy. We notice that the func¬ 
tion fi is positive definite along the attractive branch, 
whereas it changes sign at the value r = 6 on the re¬ 
pulsive branch. For positive values of b the nodal surface 
in the many-body trial wave function 07 (R) which origi¬ 


nates from the choice (22) of the Jastrow correlation term 


allows one to discriminate between the ground-state at¬ 
tractive branch and the excited-state repulsive branch. 

Furthermore, we notice that the unitary limit, corre¬ 
sponding to 6 = ±oo, is reached following the attrac¬ 
tive branch. This limit corresponds to a Jastrow term 
fi{r) oc 1/r in the range Rq < r < R and is obtained 
by approaching the resonance both from 6 > 0 and from 
& < 0 . 


fi{r) 


0 

sin[fc(r—6)] 
r 


r < b 

b < r < L/2 , 


( 21 ) 


D. Attractive and repulsive polaron branch 


similar to Eq. (20) with the only difference that the 
value of the scattering length is now given by b. 


• Square-well potential: We use different functional 
forms of fi for the repulsive and attractive branch of 
the polaron. 


Simulations are carried out using periodic boundary 
conditions and the number N of bosons in the bath is 
typically N = 64. Calculations with different numbers of 
particles up to A = 128 are also performed in order to 
check that finite-size effects are below statistical uncer¬ 
tainty. 
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Binding energy 


We determine the polaron binding energy by calculat¬ 
ing the energy difference 


^ = E{N, 1) - EoiN) , 


(24) 


where Eq{N) is the ground-state energy of the system of 
N bosons alone and E{N, 1) is the energy of the system 
of N bosons plus the impurity in the same volume V. 

Two different branches are obtained depending on the 
impurity-boson interaction potential and on the choice of 


the Jastrow term // in the trial function (19). The at 


tractive branch, which corresponds to E{N, 1) being the 
ground state of the composite system, is simulated us¬ 
ing the square-well potential (|^ and the positive definite 
function fj described in Sec. II C. Along the repulsive 
branch E{N, 1) is still the ground state of the hard-sphere 
potential (^, but it corresponds to an excited state of the 
potential El) which we calculate by imposing the nodal 


constraint given by Eq. (221 on the trial function. 


The results for both branches are shown in Fig.j^ The 
gas parameter of the bosonic bath is here na^ = 10“® 
corresponding to a dilute gas whose ground-state energy 
Eo{N) is found to be very close to the result (|^ of second 
order perturbation theory. The results reported in Fig.[^ 
are obtained both with the hard-sphere potential @ and 
with the square-well potential 0 where we used two dif¬ 
ferent values of the ratio a/Ro of the boson-boson scat¬ 
tering length to the potential range. The figure clearly 
indicates that the energies scale with the ratio a/b and 
that the details of the impurity-boson potential are ir¬ 
relevant. For the repulsive branch we find a remarkably 


good agreement with the perturbation result in Eq. (13) 


up to values of 6 ~ 25a. For larger values of the impurity- 
boson scattering length calculations using the square-well 
potential get increasingly difficult because of large fluc¬ 
tuations arising from the nodal constraint imposed on 
the Jastrow correlation term // which is no longer ad¬ 
equate to define the excited state of the polaron. It is 
worth stressing at this point that the results for the re¬ 
pulsive branch strongly depend on the choice made for 
the nodal constraint which provides the correct descrip¬ 
tion of the repulsive polaron only in the limit nb^ <C 1. 
The attractive branch can instead be followed, starting 
from small negative values of 6, down to the unitary point 
(a/6 = 0) and when approaching this point the polaron 
binding energy shows large deviations from the pertur¬ 
bation expansion ([T§. We find that at the unitary point 
1^1 ~ ‘25gn, resulting in a binding energy of the impurity 
much larger than the chemical potential of bosons in the 
bath. The results corresponding to positive values of 6 
along the attractive branch are reported in the inset of 
Fig. 0. We find that /i lies always significantly below 
the two-body binding energy Cb- One expects that, by 
increasing a/b on the positive side of the resonance, the 



FIG. 1: (color online). Polaron energy ^ as a function of the 
ratio a/b of scattering lengths for both the repulsive and the 
attractive branch. The gas parameter of the bosonic bath 
is na^ = 10~®. The symbols are the DMC results obtained 
with the following impurity-boson interaction potential: hard 
sphere (blue squares); square well with a/7?o = 5 (green 
circles); square well with a/Ro = 20 (red diamonds). The 
dashed line is the result ( |13[ ) of perturbation theory for the two 
branches. Inset: Polaron energy along the attractive branch 
on the positive side of the resonance value for the impurity- 
boson scattering length. The solid line corresponds to the 
binding energy tb in the square well potential with a/R q — 5. 


polaron binding energy eventually approaches the energy 
of the deepest cluster state (see Sec. III). 

In Fig. we show a more detailed comparison of the 
binding energy of the two polaron branches with the per¬ 


turbation expansion (13). This is carried out by sub¬ 


tracting from the values of g the mean-field contribution 


g-MF = gn-, i.e. the first term in bracket in Eq. (13), 


and by comparing g — gM f with the second order con¬ 
tribution of the perturbation expansion. We notice that 
g — gMF remains positive in the limit a/|6| <gC 1, in con¬ 
trast with the predictions of Ref. [9] and [12] where this 
quantity should turn negative for a/|6| ^ 0.07 as a con¬ 
sequence of the self-localzation of the polaron. 


• Effective mass 


The effective mass of a distinguishable particle can be 
determined in a DMC simulation by calculating its dif¬ 
fusion constant in imaginary time |341 I35j . The main 
assumption is that the energy of the system with the im¬ 
purity having momentum can be written in the form 

F(iV,l;p„)=Fo(lV)+M+^ + --- (25) 

in terms of the impurity binding energy g and effective 
mass m*. The ratio m/m* of the bare to the effective 
mass of the particle is then given by 

m (|Ar„(r)|^) 

— = lim -——- 

m* T-)-oo oDr 


( 26 ) 
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FIG. 2: (color online). Polaron energy with the mean-filed 
contribution fiMF subtracted as a function of the ratio a/b. 
The symbols are as in Fig. The dashed line is the second 
order contribution to the perturbation expansion (131 (second 
term in bracket). Notice that the results for a/b > 0 refer to 
the repulsive branch only. 


where D = It?'/ 2m is the diffusion constant of a free par¬ 
ticle and (|Ara(T)|^) = (|ra(T) — ra(0)|^) is the mean 
square displacement of the impurity in imaginary time. 
One can determine the value of m/m* from the large 
time slope of (|Ara(T)|^) as a function of the imaginary 
time T. The results are shown in Fig. [^for the attractive 
and repulsive branch of the polaron. We notice that far 
away from the resonant point a/b = 0 the increase of the 
effective mass agrees with the prediction of perturbation 
theory. On approaching the resonance, m* /m remains fi¬ 
nite reaching values ^ 2. This result is again in contrast 
with the self-localization picture of Refs. 0111] which 
predicted a large increase of the effective mass with in¬ 
creasing coupling strength. Along the attractive branch 
we calculated m* only up to the unitary point, we expect 
that following this branch on the positive side of the res¬ 
onance the value of m* should continue to increase. 


• Density profiles 


Another important output of our QMC simulations, 
useful to understand the changes induced in the bosonic 
bath by the impurity, is the pair correlation function 
gi(r) giving the probability of finding a bosonic parti¬ 
cle at a distance r from the impurity. At large distances 
gi{r) ~ 1, whereas its short-range behavior is determined 
by the impurity-boson potential Vi(r). The density pro¬ 
file of the particles of the bath surrounding the impurity 
can be calculated using the following integral of the pair 
correlation function 


n(r) 


dr'r'^gI {r') 

"" r3/3 


(27) 
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FIG. 3: (color online). Effective mass of the polaron as a 
function of the ratio a/b for the repulsive (blue circles) and 
attractive (green circles) branch. The gas parameter is given 
by na® = 10“®. The dashed line corresponds to the pertur¬ 
bation expansion (|14l). 


which approaches the bulk value n far away from the 
impurity. As a technical remark, we compute the func¬ 
tion gi(r) by carrying out both a variational and a dif¬ 
fusion Monte Carlo calculation, which provide the esti¬ 
mates gY^^ and gf^^ respectively, and by using the 
extrapolation formula ( 7 /(r) = [gf^‘"{ryY/gY^^{r) [36] . 

The ratio n{r)/n of the local to the bulk density is 
shown in Figs. 4]|6 for three values of |6|/a, both on the at¬ 
tractive and the repulsive branch of the polaron. The cal¬ 
culations are carried out at the bath density na^ = 
using the square-well potential ([^ with a/i?o = 5, and 
distances are reported in units of the healing length 
^ = l/VSTrna. In all cases n(r) exhibits a pronounced 
peak at the position of the impurity caused by the attrac¬ 
tive potential well. If 6 < 0 (attractive branch) the local 
density of particles decreases monotonously, reaching the 
bulk value n when r/^ ^ 1. Instead, if & > 0 (repulsive 
branch) the density n{r) goes through a minimum before 
reaching the bulk value. The position of the minimum 
lies in the region 0.6 < r/^ < 0.8 and decreases with in¬ 
creasing h/a. In this case, the density depletion occurring 
at a large distance from the impurity arises from the ef¬ 
fective repulsive interaction associated with the positive 
value of h. 

The average number of particles of the bath surround¬ 
ing the impurity is obtained from the integral Nb = 
Aim JY dr'r'^gj{r') and is shown in the inset of Figs, di 
[6] as a function of the distance r. We notice that the 
number Nb starts to grow faster with increasing r for 
positive values of b, consistently with the larger peak of 
n(r) at very short distances. In particular, for b/a = 10 
and 20 (see Figs. ^), Nb rapidly reaches the value of 
one particle already at distances r/^ ~ 0.1. 

We also notice that any perturbative approach based 
on the Frohlich-type Hamiltonian ( [l^ can only be mean¬ 
ingfully applied if the density perturbation induced by 
the impurity, Sn{r) = (n(r) — n), satisfies the condition 























FIG. 4: (color online). Density profile of the bath surrounding 
the impurity for hja = ±10. Inset: Integrated number of 
particles of the bath at a distance r from the impurity. 



FIG. 6: (color online). Density profile of the bath surrounding 
the impurity for hja = ±30. Inset: Integrated number of 
particles of the bath at a distance r from the impurity. 

as 


C = lim gi(r)^—{na^Y^^ 


r—>-0 


(28) 



FIG. 5: (color online). Density profile of the bath surrounding 
the impurity for hja = ±20. Inset: Integrated number of 
particles of the bath at a distance r from the impurity. 

\5n{r)\/n 1 [l^. These approaches are therefore lim¬ 

ited to values of r such that r/^ > 1 and can never de¬ 
scribe correctly the structural properties of the bath at 
short distances from the impurity. 

• Contact parameter 


where the r —>■ 0 limit should be intended in the sense 
specified above. The results for the contact parameter, 
obtained at no?' = 10“® using the SW potential with 
a/Ro = 5, are shown in Fig. for both the attractive 
and the repulsive branch. For small values of \h\/a the 
impurity-boson pair correlation function is well approxi¬ 
mated by the simple expression gi{r) = (1 — deter¬ 
mined solely by two-body physics, yielding the estimate 
C = for the contact parameter. The deriva¬ 

tive of the polaron binding energy with respect to the in¬ 
verse scattering length b should also be related to C m- 
From the behavior in the \b\/a <C 1 regime one finds 

( 29 ) 

gn d(—a/b) 

This result is shown in Fig. [^together with the contact 
extracted from the pair correlation function. Good agree¬ 
ment is found along the attractive branch, whereas the 
two estimates of C on the repulsive branch are compat¬ 
ible only in the weak-coupling limit. The disagreement 
between the contact parameter obtained from the equa¬ 
tion of state and from the pair correlation function indi¬ 
cates that our choice of the trial wave function does not 
provide a fully satisfactory description of the repulsive 
polaron in the region where b/a becomes very large. 


The impurity-boson contact parameter can be deter¬ 
mined from the behavior of the pair correlation function 
gi{r) in the range of distances r <C but still much 

larger than the typical radius of the impurity-boson inter¬ 
action. We define the dimensionless contact parameter 


E. Resonant interaction 

In this section we focus on the properties of the Bose 
polaron when the interaction between the impurity and 
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FIG. 7: (color online). Contact parameter Eq. (281 as a func¬ 
tion of the ratio ajh for the repulsive (blue circles) and at¬ 
tractive (green circles) branch. The gas parameter is given 
by na® = 10“®. The lines correspond to the determination of 
C from the attractive and repulsive branch of the equation of 
state [see Eq. (291]. 


the bath is resonant, i.e. a/h = 0. In Fig.we show the 
binding energy of the polaron calculated at resonance 
as a function of the gas parameter of the bath. The 
results show that /r scales with the energy — and 
that, once expressed in these units, it depends weakly 
on the gas parameter over many orders of magnitude. 
As no^ decreases the value of fi also decreases, reaching 
^ ~ at the very small density na^ = 3 x 10“®. 

We can not establish weather the binding energy con¬ 
tinues to decrease for even smaller densities, signaling 
the instability of the non interacting gas in the presence 
of an impurity with attractive interaction, or it reaches 
a constant value in agreement with the findings of the 
field-theoretical calculation in Ref. m- Remarkably, the 
binding energy of a Fermi polaron resonantly interact- 

t2 2/3 

ing with the bath is given by /r = —, where n 

is here the density of the Fermi sea [38], and differs ap¬ 
proximately by a factor of two compared to the results 
in Fig. for the smallest values of na^. 

The effective mass as a function of the gas parameter 
is shown in Fig. Also in this case we find a small 
variation of m*/m following a change of no? over orders 
of magnitude. The largest effective mass, m*/m ~ 1.7, 
is achieved at the smallest densities of the bath. 

In Fig. |T^ we show the density profile of the bath sur¬ 
rounding the impurity obtained using the pair correlation 


function gj and Eq. (271. The behavior is qualitatively 


similar to the one reported in Fig. [^ and corresponding 
to b/a = —30 along the attractive polaron branch. By 
decreasing the value of the bath gas parameter we find 
that the density peak around the impurity sharpens and 
the size of the deformation in units of the healing length 
shrinks. Finally, in the inset of Fig. we show the 
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FIG. 8: Polaron binding energy at unitarity (a/6 = 0) in units 

t2 2/3 

of — as a function of the gas parameter of the bath. 



FIG. 9: Effective mass of the polaron at unitarity (a/6 = 0) 
as a function of the gas parameter of the bath. 

value of the contact parameter C, determined from the 
short-range behavior of the pair correlation function, for 
different values of Also for this quantity we 

observe a weak dependence on the value of the bath gas 
parameter. 


III. FEW-BODY PHYSICS 


In this Section we consider the problem of the existence 
of bound states in vacuum consisting of the impurity and 
a number N of bosons. Of course, such bound states can 
only occur in the case of the square-well model for the 
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FIG. 10: Density profile of the bath surrounding the impu¬ 
rity at unitarity (o/6 = 0) for three values of the gas parame¬ 
ter. Inset: Contact parameter C at unitarity as a function of 

(na3)i/3. 



FIG. 11: Energy of clusters with A^-l-l particles at the unitary 
point {a/b = 0) as a function of the number N of bosons. 

t 2 

Energies are in units of value at = 1 refers to 

the two-body binding energy. Inset: Same as main figure at 
a/h — 0.1. 


impurity-boson potential. This potential supports a two- 
body molecular state having energy e;,, given by Eq. (§, 
for all positive values of the inter-species scattering length 
b. 

The search for the ground state of clusters with IV -|-1 
particles is carried out using the DMC method based on 
the following trial wave function 


N 


tprCR) = expi-^Rn) //(Ca) /b(Gj) . (30) 

i=l i<j 


The above wave function differs from the one of Eq. (191, 
used in simulations of homogeneous configurations, by 


the exponential term which depends on the hyper-radius 
of the cluster 


Rh = 


\ 


N 


- TCm)^ + “ ^CmY 


(31) 


2=1 


where vcm = (la w 
the center of mass. The Jastrow correlation terms in 


2 ^ 2=1 


is the coordinate of 


Eq. (30) are similar to the ones of Eqs. (20) and (23), 


respectively for the boson-boson fs and the impurity- 
boson fj function. Since periodic boundary conditions 


are absent here, the length scale L/2 in Eq. (20) is re¬ 


placed by the large distance Rl = 300a. Moreover, the 
boundary condition on the derivative of fj is relaxed with 
the choice fi{r) = B + holding for r > R, with 

the constants B and C determined in the same way as 


in Eq. (23). Free parameters that are optimized using a 


variational procedure are the matching point R and the 
coefficients a and /?. In particular, the latter fixes the 
size of the cluster in terms of its hyper-radius. 

Calculations are performed in the reference frame 
where ycm — 0, in order to eliminate the contribution 
from the center of mass motion. Furthermore, we con¬ 
sider only the resonant point where a/b = 0 and Cb = 0, 
and the point a/b = 0.1 on the positive side of the reso¬ 
nance where Cb < 0. In Fig.[^we show the results for the 
ground-state energy of the cluster with iV -|- 1 particles 
as a function of the number N of bosons. At unitarity 
the two-body binding energy, corresponding to = 1 in 
Fig.[TTJ is identically zero, whereas the three-body Efi¬ 
mov state {N = 2) is found to feature an extremely shal¬ 
low ground-state energy: IE 2 + 1 I ^ 10~'^ 2ma3 • This result 
is consistent with the prediction E 2+1 = — 

mrQ 

for the lowest Efimov state in terms of the three-body 
length To and the Efimov parameter sg |39j . In the 
case of equal masses for the impurity and the bosons, 
the value of Sg is very small, sg = 0.4137, resulting in 
E 2+1 ^ -10“°-^ which is of the same order as our esti- 

mrQ 

mate if rg ~ a. For increasing N the cluster ground-state 
energy decreases markedly up to fV = 4, while clusters 
with N = 6 are undoubtedly unbound. These findings 
are compatible with the results at a/b = 0.1 (see inset 
of Fig. 11), where the binding energy appears not to de¬ 


crease further already for N > 3. 

It is important to stress that the largest size N of 
bound clusters, as well as the precise value of their 
ground-state energies, depend on the details of the inter¬ 
boson and impurity-boson interactions. However, we be¬ 
lieve that the qualitative behavior emerging from our 
simulations should hold for any short-range interaction 
with scattering length a and b, respectively. In particu¬ 
lar, we notice that the polaron binding energy shown in 
Fig.[T] is more than a factor 10^ larger than the deep¬ 
est cluster state at unitarity and remains larger also at 
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ajh = 0.1. A possible reason for the irrelevance of clus¬ 
ter states at unitarity is the feature of equal masses for 
the impurity and the bosons, which makes Efimov states 
extremely shallow. 


IV. MANY IMPURITIES 


Let us now analyze the case of a small concentration 
of impurities immersed in a BEG at T = 0. For this 
problem the statistics of the impurities is important and 
in the present study we consider only impurities which 
obey Bose statistics. A binary Bose-Fermi mixture with 
a small concentration of bosons in a Fermi sea and fea¬ 
turing resonant Bose-Fermi interactions has been inves¬ 
tigated using QMC methods in Ref. [40| . 

A collection of M impurities immersed in a gas of N 
particles is described by the Hamiltonian 

^2 N 

''<>('■«) (32) 

i—1 i<j 

^2 M N M 

- w E ^« + E >'"<''<.<!) + E E ‘''f’'*") ■ 

a —1 q;<^ 2—1 a—1 


The impurity-impurity potential is modeled by the same 
hard-sphere interaction, including the same scattering 
length a, which characterizes the coupling between the 
bosons of the bath: U//(r) = Vsir). We also assume 
that the masses of the two types of particles are the same 
(mi = ttib), and the impurity-boson interaction Vi(r) is 
as described in Sec. II-A. 

The perturbation treatment of a binary mixture of 
Bose condensates at T = 0 has been carried out in 
Ref. [41] using an extension of the standard Bogoliubov 
approach. The result for the ground-state energy in the 
low concentration limit, x = M/N ^ 1, is obtained as 
follows 


Eq = Eb 


Ngn 


1-b 


b 

a 

64 


32 /—6^ 
,2 


3v^ a2 ) 2 


(33) 


up to quadratic contributions in the ratio b/a of scatter¬ 
ing lengths and in the impurity concentration. Here Eb 
is the energy ^ of the bath without impurities. Fur¬ 
thermore, one should notice that the term linear in the 
concentration coincides with the single polaron energy of 
Eq. (131, whereas the term proportional to describes 
the repulsive interaction between polarons. 

In QMC simulations we calculate the ground-state en¬ 
ergy E{N, M) of the mixture of N bosons plus M impu¬ 
rities making use of the following trial wave function 


N M 

^t(R.) = W fB(r,j) n feirap) H IT /KUa) , (34) 

i<j a<l3 2=1q;=1 



FIG. 12: Energy shift between the system with a given con¬ 
centration X of impurities and the bath without impurities as 
a function of the ratio of scattering lengths along the repulsive 
branch. The density of the bath is na^ — 10~® and results 
are shown for three different concentrations. Solid lines refer 
to the perturbation result Eq — Eb from Eq. (|33[). 


where the same Jastrow factor /s of Eq. (21) accounts 
for the repulsive correlations between the particles of 
the bath and between the impurities while the impurity- 
boson term // is as described in Sec. H-C. We notice that 
the above wave function is symmetric under the exchange 
of the impurity coordinates fulfilling Bose statistics. 

We perform QMC calculations using N = 64 particles 
in the bath and a varying number M < 15 of impurities 
in a cubic box with periodic boundary conditions, aiming 
to simulate a homogeneous binary mixture characterized 
by a small concentration, x < 0.23, of the minority com¬ 
ponent. In particular, we calculate the shift between the 
ground-state energy of the mixture and of the bath with¬ 
out impurities: 


AF;(V, M) = E{N, M) - EoiN) . 


(35) 


The results for different values of the ratio b/a along the 
repulsive branch are reported in Fig. for the bath den¬ 
sity na^ = 10“®. The energy shift AE is in remarkable 


agreement with the perturbation result (331 for all values 


of x up to b/a ~ 8. Only at the largest values of b/a and 


of the concentration signihcant deviations from Eq. (33) 


are visible. In particular, the interaction between po¬ 
larons appears to be overestimated by the term propor¬ 
tional to x^ in Eq. (33). In Fig. 13 we show the results 


of AE as a function of the concentration x for b/a = ±5 
on both the repulsive and the attractive branch. Also 


in this case, the comparison with Eq. (33) shows a very 
good agreement. 


The above results validate the expression (33) for the 


equation of state of a binary Bose-Bose mixture in the 
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condition Vi + V 2 = V. The stability of this state re¬ 


quires that n = m and that the energy cost to add one 


impurity to each of the two uniform ph ases is positive: 
^ > 0- From Eqs. ^ and (13), giving respec¬ 

tively the energy of the bath without impurities and the 


excess energy of a single polaron, we get that the phase 
separated state is stable if 




a 



if 6 > 0 (38) 

if 6 < 0 . (39) 


One can easily show that, if 6 > 0, the energy (37) of 
the phase separated state lies below the energy (33) of 
the homogeneous binary mixture for any value ^ > 1. 
On the contrary, no gas-like phase appears to be stable 
outside the regions of Eqs. (38) and (36) when b <0. 


FIG. 13: Energy shift AE{M, N) as a function of the concen¬ 
tration of impurities for the two values ^ = ±5 of the ratio 
of scattering lengths. The density of the bath is as in Fig. |12[ 
Dashed lines refer to the perturbation result Eq — Eb from 

Eq. pi. 


regime of parameters: \b\/a < 5 and x < 0.2. Such a 
validation is important in order to establish the stabil¬ 
ity conditions and the phase diagram of the mixture in a 
quantitatively reliable way. From the analysis of the com¬ 
pressibility matrix ^ ^ ^ where ni = N/V 

and 712 = M/V, one finds that the homogeneous mix¬ 
ture is stable if the ratio of scattering length satisfies the 
condition 


- 1 - 



'na^ < - < 1 -I- 



(36) 


This relation holds in the limit a; <C 1 and, compared to 
the mean-field result given by |&|/a < 1 [H], includes the 
leading order correction in the small parameter V na^. 
We notice that the results reported in Figs. Ill and Ill 


lie outside the stability range (36) of the homogeneous 
binary mixture. The spinodal instability arising from 
a vanishing compressibility is associated to systems ap¬ 
proaching the thermodynamic limit and is usually pre¬ 
vented in simulations of finite-size systems. 

Another possible state of the binary mixture corre¬ 
sponds to a complete phase separation between the N 
bosons and the M impurities. In this case the ground- 
state energy can be written as 


E{ 


\ TV 

n,m) = n — 


-b V2 



(37) 


in terms of the densities n = N/Vi and m = M/V 2 of 
the two species and their relative volumes fulfilling the 


V. CONCLUSIONS 

We investigated the properties of an impurity im¬ 
mersed in a Bose-Einstein condensate at T = 0. This 
Bose polaron study has been carried out using a fully 
microscopic approach and QMC simulation methods in 
analogy with previous investigations of the more thor¬ 
oughly expounded, both theoretically and experimen¬ 
tally, Fermi polaron problem. The main results concern 
the binding energy of the impurity and its effective mass 
along the attractive and repulsive polaron branch, ex¬ 
plored by changing the scattering length of the impurity- 
boson interaction potential. These results are expected 
to be universal, for a given boson-boson and impurity- 
boson scattering length and for equal masses of the two 
components. At the resonant point of the impurity-boson 
interaction the polaron binding energy scales with the 
equivalent of the Fermi energy of the bath in analogy 
with the behavior found for Fermi polarons. 

The measurement of the polaron binding energy should 
be accessible in experiments using radio-frequency spec¬ 
troscopy while its effective mass affects the disper¬ 
sion of collective excitations in highly imbalanced two- 
component mixtures. We hope that such experiments, 
which were so successful in the investigation of the prop¬ 
erties of the Fermi polaron, will be carried out also for its 
Bose counterpart providing a deeper knowledge of this 
clean and simply stated, but highly non trivial many- 
body problem. 
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